﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Diplomova_prace
{
    class Linear_Solutions_alg
    {
        LinearniRovnice_ret ret = new LinearniRovnice_ret();
        LinearniRovniceSlozky_ret ret1 = new LinearniRovniceSlozky_ret();
        double scalarDet = 1; // pro vypocet je potreba porovnani s DP pred 9.3.2014 //Gaussova eliminace 
        public double[] backTracking(double[,] matrix, double[,] vectorB)
        {
            int rows = matrix.GetLength(0);
            int columns = matrix.GetLength(1)-1;
            double[] vectorX = new double[rows];
            for (int i = rows; i > -1; i--)
            {
                for (int j = columns; j > -1; j--)
                {
                    if (j > i)
                        vectorB[i,0] -= matrix[i, j] * vectorX[j];
                    else if (j == i)
                        vectorX[i] = vectorB[i,0] / matrix[i, j];
                    else
                        break;
                }
            }
            return vectorX;
        }
        public double[,] pivoting(double[,] matrix, int column,double[,] vectorB, double[,] matrixP = null)
        {
            double pivot = matrix[column, column];
            int change = column;
            double[] tmp;
            double temp;
            for (int i = column; i < matrix.GetLength(0); i++)
            {
                if (Math.Abs(matrix[i, column]) > Math.Abs(pivot))
                {
                    pivot = matrix[i, column];
                    change = i;
                }
            }
            if (change != column) 
            {
                tmp = new double[matrix.GetLength(1)];
                for (int i = 0; i < matrix.GetLength(1); i++)
                {
                    tmp[i] = matrix[column, i];
                    matrix[column, i] = matrix[change, i];
                    matrix[change, i] = tmp[i];
                }
               
            }
            if (vectorB != null)
            {
                temp = vectorB[change, 0];
                vectorB[change, 0] = vectorB[column, 0];
                vectorB[column, 0] = temp;
            }
            if (matrixP != null)
            {
                tmp = new double[matrix.GetLength(1)];
                for (int i = 0; i < matrixP.GetLength(1); i++)
                {
                    tmp[i] = matrixP[column, i];
                    matrixP[column, i] = matrixP[change, i];
                    matrixP[change, i] = tmp[i];
                }
            }
            return matrix;

        }
        public double[,] Ly(double[,] matrixL, double[,] vectorB,double[,] matrixP=null)
        {
            int rows = matrixL.GetLength(0);
            int columns = matrixL.GetLength(1);
            double[,] vectors = new double[vectorB.GetLength(0),1];
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < columns; j++)
                {
                    if (j < i)
                        vectorB[i,0] -= matrixL[i, j] * vectors[j,0];
                    else if (j == i)
                        vectors[i,0] = vectorB[i,0] / matrixL[i, j];
                    else
                        break;
                }
            }
            return  vectors;
        }
        public double[,] Ux(double[,] matrixU, double[,] vectors)
        {
            int rows = matrixU.GetLength(0);
            int columns = matrixU.GetLength(1);
            double[,] vectorX = new double[vectors.Length,1];
            for (int i = rows-1; i > -1; i--)
            {
                for (int j = columns-1; j > -1; j--)
                {
                    if (j > i)
                        vectors[i,0] -= matrixU[i, j] * vectorX[j,0];
                    else if (j == i)
                        vectorX[i,0] = vectors[i,0] / matrixU[i, j];
                    else
                        break;
                }
            }
            return vectorX;
        }
        
        public Tuple<List<double[,]>,List<double[,]>,List<double>> JacobihoMethodMatrix(double[,] matrix, double[,] vectorB,double[,] vectorX, int countOfIteration,double accuracy)
        {
            List<double> ListAccuracy = new List<double>();
            List<double[,]> vectorsX = new List<double[,]>() {vectorX };
            double[,] matrixE = new double[matrix.GetLength(0), matrix.GetLength(1)];
            double[,] matrixDInverse = new double[matrix.GetLength(0), matrix.GetLength(1)];
            double[,] matrixF = new double[matrix.GetLength(0), matrix.GetLength(1)];
              
            for (int i = 0; i < matrix.GetLength(0); i++)
            {
                for (int j = 0; j < matrix.GetLength(1); j++)
                {
                    if (i == j)
                    {
                        matrixDInverse[i,j] = 1 / matrix[i, j];
                    }
                    else if ( i < j)
                    {
                        matrixF[i,j] =  matrix[i, j];
                    }
                    else
                    {
                        matrixE[i, j] = matrix[i, j];
                    }
                }
            }
            double[,] tmpRightSide = MulOFMatrix(matrixDInverse, vectorB);
            double[,] EplusF = SumOfMatrix(matrixE, matrixF);
            double[,] MultiplyWithConstanto = MultiplyWithConstant(matrixDInverse, -1);
            double[,] tmp = MulOFMatrix(MultiplyWithConstanto,EplusF);
            double[,] temp;
            for (int i = 0; i < countOfIteration; i++)
            {
                temp = Copy(vectorX);
                
                vectorX = SumOfMatrix(MulOFMatrix(tmp, vectorX),tmpRightSide);
                vectorsX.Add(vectorX);
                ListAccuracy.Add(EndingConditionOfIteration(vectorX, temp));
                if (Math.Abs(ListAccuracy.Last()) < accuracy)
                {
                    return Tuple.Create(new List<double[,]>() { matrixDInverse, matrixE, matrixF, vectorB, vectorsX[0], EplusF, tmp, tmpRightSide }, vectorsX, ListAccuracy);
                }
            }
            return Tuple.Create(new List<double[,]>() { matrixDInverse, matrixE, matrixF, vectorB, vectorsX[0], EplusF, tmp, tmpRightSide }, vectorsX,ListAccuracy);
        }
        public Tuple<List<double[,]>, List<double[,]>, List<double>> GaussSeidelMethodMatrix(double[,] matrix, double[,] vectorB, double[,] vectorX, int countOfIteration, double accuracy)
        {
            double[,] matrixE = new double[matrix.GetLength(0), matrix.GetLength(1)];
            double[,] matrixD = new double[matrix.GetLength(0), matrix.GetLength(1)];
            double[,] matrixF = new double[matrix.GetLength(0), matrix.GetLength(1)];
            List<double> Listaccuracy = new List<double>();
            List<double[,]> vectorsX = new List<double[,]>() { vectorX };

            for (int i = 0; i < matrix.GetLength(0); i++)
            {
                for (int j = 0; j < matrix.GetLength(1); j++)
                {
                    if (i == j)
                    {
                        matrixD[i, j] = matrix[i, j];
                    }
                    else if (i < j)
                    {
                        matrixF[i, j] = matrix[i, j];
                    }
                    else
                    {
                        matrixE[i, j] = matrix[i, j];
                    }
                }
            }
            double[,] InversionDplusE = InversionOfMatrix(SumOfMatrix(matrixD, matrixE));
            double[,] tmpRightSide = MulOFMatrix(InversionDplusE, vectorB);
            double[,] tmp = MulOFMatrix(MultiplyWithConstant(InversionDplusE, -1),  matrixF);
            double[,] temp;
            for (int i = 0; i < countOfIteration; i++)
            {
                temp = Copy(vectorX);
                vectorX = SumOfMatrix(MulOFMatrix(tmp, vectorX), tmpRightSide);
                vectorsX.Add(vectorX);
                Listaccuracy.Add(EndingConditionOfIteration(vectorX, temp));
                if (Math.Abs(Listaccuracy.Last()) < accuracy)
                {
                    return Tuple.Create(new List<double[,]>() { matrixD, matrixE, matrixF, vectorB, vectorsX[0], SumOfMatrix(matrixD, matrixE), InversionDplusE, tmp, tmpRightSide }, vectorsX, Listaccuracy);
                }
            }
            return Tuple.Create(new List<double[,]>() { matrixD, matrixE, matrixF, vectorB, vectorsX[0], SumOfMatrix(matrixD, matrixE), InversionDplusE, tmp, tmpRightSide }, vectorsX ,Listaccuracy);
        }
        public Tuple<List<double[,]>, List<double[,]>,List<double>> SORMethodMatrix(double[,] matrix, double[,] vectorB, double factor, double[,] vectorX,int countOfIteration, double accuracy)
        {
            double[,] matrixE = new double[matrix.GetLength(0), matrix.GetLength(1)];
            double[,] matrixD = new double[matrix.GetLength(0), matrix.GetLength(1)];
            double[,] matrixF = new double[matrix.GetLength(0), matrix.GetLength(1)];
            List<double> Listaccuracy = new List<double>();
            List<double[,]> vectorsX = new List<double[,]>() { vectorX };
       
            for (int i = 0; i < matrix.GetLength(0); i++)
            {
                for (int j = 0; j < matrix.GetLength(1); j++)
                {
                    if (i == j)
                    {
                        matrixD[i, j] = matrix[i, j];
                    }
                    else if (i < j)
                    {
                        matrixF[i, j] = matrix[i, j];
                    }
                    else
                    {
                        matrixE[i, j] = matrix[i, j];
                    }
                }
            }
            double[,] wF = MultiplyWithConstant(matrixF, factor);
            double[,] wE = MultiplyWithConstant(matrixE, factor);
            double[,] DpluswE = SumOfMatrix(matrixD, wE);
            double[,] InversionDpluswE = InversionOfMatrix(DpluswE);
            double[,] JednaMinuswD = MultiplyWithConstant(matrixD, (1 - factor));
            double[,] tmpLeftSide = DiffOfMatrix(JednaMinuswD, wF);
            double[,] wB = MultiplyWithConstant(vectorB, factor);
            double[,] tmpRightSide = MulOFMatrix(InversionDpluswE, wB);
            double[,] tmp = MulOFMatrix(InversionDpluswE, tmpLeftSide);
            for (int i = 0; i < countOfIteration; i++)
            {
                double[,] temp = Copy(vectorX);
                vectorX = SumOfMatrix(MulOFMatrix(tmp, vectorX), tmpRightSide);
                vectorsX.Add(vectorX);
                Listaccuracy.Add(EndingConditionOfIteration(vectorX, temp));
                if (Math.Abs(Listaccuracy.Last()) < accuracy)
                {
                    return Tuple.Create(new List<double[,]>() { matrixD, matrixE, matrixF, vectorB, vectorsX[0], wE, wF, wB,DpluswE, InversionDpluswE, JednaMinuswD, tmp, tmpRightSide, }, vectorsX, Listaccuracy);
                }
            }
                return Tuple.Create(new List<double[,]>() { matrixD, matrixE, matrixF, vectorB, vectorsX[0], wE, wF, wB,DpluswE, InversionDpluswE, JednaMinuswD, tmp, tmpRightSide, }, vectorsX, Listaccuracy );
        }
        public LinearniRovniceSlozky_ret SORMethodFolder(double[,] matrix, double[,] vectorB, double[,] vectorX, double factor, int countOfIteration, double accuracy)
        {
            List<double> Listaccuracy = new List<double>();
            List<double[,]> res = new List<double[,]>() { Copy(vectorX) };
            double[,] temp = new double[matrix.GetLength(0), 1];
            double Sum;
            for (int i = 0; i < countOfIteration; i++)
            {
                temp = Copy(vectorX);
                for (int j = 0; j < matrix.GetLength(0); j++)
                {
                    Sum = 0;
                    for (int k = 0; k < matrix.GetLength(1); k++)
                    {
                        if (j != k)
                            Sum += matrix[j, k] * vectorX[k, 0];
                    }
                    Sum -= vectorB[j, 0];
                    vectorX[j, 0] = (1 - factor) * vectorX[j, 0] - ((factor / matrix[j, j]) * Sum);
                }

                Listaccuracy.Add(EndingConditionOfIteration(vectorX, temp));
                res.Add(Copy(vectorX));
                if (Math.Abs(Listaccuracy.Last()) < accuracy)
                {
                    ret1.vectorBProperty = vectorB;
                    ret1.vectorXProperty = vectorX;
                    ret1.matrixProperty = matrix;
                    ret1.vectorsProperty = res;
                    ret1.factorProperty = factor;
                    ret1.vectorsAccuracyProperty = Listaccuracy;
                    return ret1;
                }
            }
            ret1.vectorBProperty = vectorB;
            ret1.vectorXProperty = vectorX;
            ret1.matrixProperty = matrix;
            ret1.vectorsProperty = res;
            ret1.factorProperty = factor;
            ret1.vectorsAccuracyProperty = Listaccuracy;
            return ret1;
        }
        public LinearniRovniceSlozky_ret GaussSeidelMethodFolder(double[,] matrix, double[,] vectorB, double[,] vectorX, int countOfIteration, double accuracy)
        {
            List<double> Listaccuracy = new List<double>();
            List<double[,]> res = new List<double[,]>() { Copy(vectorX) };
            double[,] temp = new double[matrix.GetLength(0), 1];
            for (int i = 0; i < countOfIteration; i++)
            {
                temp = Copy(vectorX);
                for (int j = 0; j < matrix.GetLength(0); j++)
                {
                    double Sum = 0;
                    for (int k = 0; k < matrix.GetLength(1); k++)
                    {
                        if (j != k)
                            Sum += matrix[j, k] * vectorX[k, 0];
                    }
                    vectorX[j, 0] = (-1 / matrix[j, j]) * Sum + (vectorB[j, 0] / matrix[j, j]);
                }
                Listaccuracy.Add(EndingConditionOfIteration(vectorX, temp));
                res.Add(Copy(vectorX));
                if (Math.Abs(Listaccuracy.Last()) < accuracy)
                {

                    ret1.vectorBProperty = vectorB;
                    ret1.vectorXProperty = vectorX;
                    ret1.matrixProperty = matrix;
                    ret1.vectorsProperty = res;
                    ret1.vectorsAccuracyProperty = Listaccuracy;
                    return ret1;
                }
            }
            ret1.vectorBProperty = vectorB;
            ret1.vectorXProperty = vectorX;
            ret1.matrixProperty = matrix;
            ret1.vectorsProperty = res;
            ret1.vectorsAccuracyProperty = Listaccuracy;
            return ret1;

        }
        public LinearniRovniceSlozky_ret JacobihoMethodFolder(double[,] matrix, double[,] vectorB, double[,] vectorX, int countOfIteration, double accuracy)
        {
            List<double[,]> res = new List<double[,]>() { Copy(vectorX) };
            List<double> Listaccuracy = new List<double>();
            double[,] temp = new double[matrix.GetLength(0), 1];
            double Sum;
            for (int i = 0; i < countOfIteration; i++)
            {
                temp = Copy(vectorX);
                for (int j = 0; j < matrix.GetLength(0); j++)
                {

                    Sum = 0;
                    for (int k = 0; k < matrix.GetLength(1); k++)
                    {
                        if (j != k)
                            Sum += matrix[j, k] * temp[k, 0];
                    }
                    vectorX[j, 0] = (-1 / matrix[j, j]) * Sum + (vectorB[j, 0] / matrix[j, j]);
                }
                Listaccuracy.Add(EndingConditionOfIteration(vectorX, temp));
                res.Add(Copy(vectorX));
                if (Math.Abs(Listaccuracy.Last()) < accuracy)
                {
                    ret1.vectorBProperty = vectorB;
                    ret1.vectorXProperty = vectorX;
                    ret1.matrixProperty = matrix;
                    ret1.vectorsProperty = res;
                    ret1.vectorsAccuracyProperty = Listaccuracy;
                    return ret1;
                }
            }
            ret1.vectorBProperty = vectorB;
            ret1.vectorXProperty = vectorX;
            ret1.matrixProperty = matrix;
            ret1.vectorsProperty = res;
            ret1.vectorsAccuracyProperty = Listaccuracy;
            return ret1;
        }
        
        public LinearniRovnice_ret CholeskyDecomposition(double[,] matrix, double[,] vectorB)
        {
            double[,] matrixL = new double[matrix.GetLength(0), matrix.GetLength(1)];
            List<double[,]> diagonal = new List<double[,]>();
            double Sum = 0;
            for (int i = 0; i < matrix.GetLength(1); i++)
            {
                for (int j = 0; j < matrix.GetLength(0); j++)
                {
                    if (j == i)
                    {
                        Sum = 0;
                        for (int k = 0; k < i; k++)
                        {
                            Sum += Math.Pow(matrixL[j, k], 2);
                        }
                        matrixL[j, j] = Math.Sqrt(matrix[j, j] - Sum);
                    }
                    if (j >= i)
                    {
                        Sum = 0;
                        for (int k = 0; k < i; k++)
                        {
                            Sum += matrixL[j, k] * matrixL[i, k];
                        }
                        matrixL[j, i] = (matrix[j, i] - Sum) / matrixL[i, i];
                    }
                }
            }
            double[,] matrixLTran = Transposition(matrixL);
            double[,] vB = Copy(vectorB);
            double[,] vectors = Ly(matrixL, vectorB);
            double[,] vY = Copy(vectors);
            double[,] vectorX = Ux(matrixLTran, vectors);

            ret.matrixLProperty = matrixL;
            ret.matrixUProperty = matrixLTran;
            ret.matrixProperty = matrix;
            ret.vectorBProperty = vB;
            ret.vectorYProperty = vY;
            ret.vectorXProperty = vectorX;
            return ret;
        }
        public LinearniRovnice_ret LUDecomposition(double[,] matrix, double[,] vectorB)
        {
            int rows = matrix.GetLength(0);
            int columns = matrix.GetLength(1);
            double[,] matrixL = new double[rows, columns];
            double[,] matrixU = new double[rows, columns];
            double[,] matrixP = new double[rows, columns];
            for (int i = 0; i < rows; i++)
            {
                matrixL[i, i] = 1;
                matrixP[i, i] = 1;
            }
            for (int i = 0; i < rows; i++)
            {
                pivoting(matrix, i, null, matrixP);
            }

            double sumU = 0;
            double sumL = 0;
            for (int i = 0; i < columns; i++)
            {
                for (int j = 0; j < i + 1; j++)
                {
                    sumU = 0;
                    for (int k = 0; k < j; k++)
                    {
                        sumU += matrixL[j, k] * matrixU[k, i];
                    }
                    matrixU[j, i] = matrix[j, i] - sumU;
                }

                for (int j = i + 1; j < columns; j++)
                {
                    sumL = 0;
                    for (int k = 0; k < i; k++)
                    {
                        sumL += matrixL[j, k] * matrixU[k, i];
                    }
                    matrixL[j, i] = (1 / matrixU[i, i]) * (matrix[j, i] - sumL);
                }
            }
            ret.vectorBProperty = Copy(MulOFMatrix(matrixP, vectorB));
            ret.matrixPProperty = Transposition(matrixP);
            vectorB = Copy(MulOFMatrix(matrixP, vectorB));
            double[,] vectors = Ly(matrixL, vectorB, matrixP);
            ret.vectorYProperty = Copy(vectors);
            ret.vectorXProperty = Ux(matrixU, vectors);
            ret.matrixProperty = matrix;
            ret.matrixLProperty = matrixL;
            ret.matrixUProperty = matrixU;
            return ret;
        }
        public Tuple<List<double[,]>, List<double[,]>,List<double>, double[,]> GaussElimination(double[,] matrixA, double[,] vectorB = null)
        {
            List<double[,]> resM = new List<double[,]>();
            List<double[,]> resV = new List<double[,]>();
            double[,] matrix = Copy(matrixA);
            double[,] vector;
            if (vectorB != null)
            {
                vector = Copy(vectorB);
            }
            else
            {
                vector = new double[matrix.GetLength(0), 1];
            }

            resM.Add(Copy(matrix));
            resV.Add(Copy(vector));
            int PivotingColumn = 0;
            resM.Add(Copy(pivoting(matrix, PivotingColumn, vector, null)));
            resV.Add(Copy(vector));
            double[,] lastM = Copy(matrix);
            double[,] lastV = Copy(vector);
            int CurrentEditRow = 1;
            double[,] workingMatrix = new double[matrix.GetLength(0), matrix.GetLength(1)];
            double[,] workingVector = new double[matrix.GetLength(0), 1];
            double ratio = double.NaN;
            int count = CountPassing(matrix.GetLength(0) - 1);
            int increment = matrix.GetLength(0) - 1;
            int stop = increment - 1;
            int ZerosBeforeDiagonal = 1;
            int offset = 1;
            bool breakChecking = true;

            double[,] tmp = new double[matrix.GetLength(0), matrix.GetLength(1)];
            List<double> ratios = new List<double>();
            for (int k = 0; k < count; k++)
            {
                for (int i = 0; i < matrix.GetLength(0); i++)
                {
                    for (int j = 0; j < matrix.GetLength(1); j++)
                    {
                        if (i != CurrentEditRow)
                        {
                            workingMatrix[i, j] = lastM[i, j];
                            workingVector[i, 0] = lastV[i, 0];
                        }
                        else
                        {
                            if (j < ZerosBeforeDiagonal)
                            {
                                if (j + 1 == ZerosBeforeDiagonal)
                                {
                                    ratio = lastM[i, j] / lastM[j, j];
                                    ratios.Add(ratio);
                                }
                            }
                            else
                            {
                                workingMatrix[i, j] = lastM[i, j] - ratio * lastM[i - offset, j];
                                workingVector[i, 0] = lastV[i, 0] - ratio * lastV[i - offset, 0];
                            }
                        }
                    }
                }
                for (int i = 0; i < matrix.GetLength(0); i++)
                {
                    for (int j = 0; j < matrix.GetLength(1); j++)
                    {
                        if (i > j && workingMatrix[i, j] != 0)
                        {
                            breakChecking = false;
                            break;
                        }
                    }
                    if (!breakChecking)
                        break;
                }
                if (breakChecking)
                {
                    resM.Add(Copy(workingMatrix));
                    resV.Add(Copy(workingVector));
                    return Tuple.Create(resM, resV, ratios, Ux(Copy(resM.Last()), Copy(resV.Last())));
                }
                else
                    breakChecking = true;
                resM.Add(Copy(workingMatrix));
                resV.Add(Copy(workingVector));
                offset++;
                CurrentEditRow++;
                if (k == stop)
                {
                    increment--;
                    stop = k + increment;
                    ZerosBeforeDiagonal++;
                    CurrentEditRow = ZerosBeforeDiagonal;
                    offset = 1;
                    PivotingColumn++;
                    resM.Add(Copy(pivoting(workingMatrix, PivotingColumn, workingVector, null)));
                    resV.Add(Copy(workingVector));
                    
                }
                lastM = Copy(workingMatrix);
                lastV = Copy(workingVector);
                workingVector = new double[matrix.GetLength(0), 1];
                workingMatrix = new double[matrix.GetLength(0), matrix.GetLength(1)];
            }
            
            return Tuple.Create(resM, resV, ratios, Ux(Copy(resM.Last()), Copy(resV.Last())));
        }

        public int CountPassing(int num)
        {
            int res = 0;
            for (int i = num; i > 0; i--)
                res += i;
            return res;
        }
        public Tuple<double[,],double[,]> GaussEliminationSupport(double[,] matrixA, double[,] vectorB=null)
        {
            double[,] matrix = Copy(matrixA);
            double[,] vector;
            if (vectorB != null)
            {
                vector = Copy(vectorB);
            }
            else
            {
                vector = new double[matrix.GetLength(0), 1];
            }
            int CurrentEditRow = 1;
            double ratio = double.NaN;
            int count = CountPassing(matrix.GetLength(0) - 1);
            int increment = matrix.GetLength(0) - 1;
            int stop = increment - 1;
            int ZerosBeforeDiagonal = 1;
            int offset = 1;
            for (int k = 0; k < count; k++)
            {
                for (int i = 0; i < matrix.GetLength(0); i++)
                {
                    for (int j = 0; j < matrix.GetLength(1); j++)
                    {
                        if (i != CurrentEditRow)
                        {
                            break;
                        }
                        else
                        {
                            if (j < ZerosBeforeDiagonal)
                            {
                                if (j + 1 == ZerosBeforeDiagonal)
                                {
                                    ratio = matrix[i, j] / matrix[j, j];
                                }
                                matrix[i, j] = 0;
                            }
                            else
                            {
                                matrix[i, j] = matrix[i, j] - ratio * matrix[i - offset, j];
                                if (j + 1 == matrix.GetLength(1))
                                    vector[i, 0] = vector[i, 0] - ratio * vector[i - offset, 0];
                            }
                        }
                    }
                }
                offset++;
                CurrentEditRow++;
                if (k == stop)
                {
                    increment--;
                    stop = k + increment;
                    ZerosBeforeDiagonal++;
                    CurrentEditRow = ZerosBeforeDiagonal;
                    offset = 1;
                    for (int i = 0; i < matrix.GetLength(0); i++)
                    {
                        if (matrix[i, i] == 0)
                            matrix = pivoting(matrix, i, vector, null);
                    }
                }
            }
            return Tuple.Create(matrix, vector);
        }
        public double[,] Copy(double[,] matrixA)
           {
               double[,] matrixB = new double[matrixA.GetLength(0), matrixA.GetLength(1)];
               for (int i = 0; i < matrixA.GetLength(0); i++)
               {
                   for (int j = 0; j < matrixA.GetLength(1); j++)
                   {
                       matrixB[i, j] = matrixA[i, j];
                   }
               }
               return matrixB;
           }
        public double[,] DiffOfMatrix(double[,] matrixA, double[,] matrixB)
        {
             double[,] matrix = new double[matrixA.GetLength(0),matrixB.GetLength(1)];

                for (int i = 0; i < matrixA.GetLength(0); i++)
                {
                   
                    for (int j = 0; j < matrixA.GetLength(1); j++)
                    {
                        matrix[i, j] = matrixA[i, j] - matrixB[i, j];
                    }
                }
            return matrix;
        }
        public double[,] SumOfMatrix(double[,] matrixA, double[,] matrixB)
        {
            double[,] matrix = new double[matrixA.GetLength(0), matrixB.GetLength(1)];

            for (int i = 0; i < matrixA.GetLength(0); i++)
            {

                for (int j = 0; j < matrixA.GetLength(1); j++)
                {
                    matrix[i, j] = matrixA[i, j] + matrixB[i, j];
                }
            }
            return matrix;
        }
        public double[,] MultiplyWithConstant(double[,] matrixA, double konstanta)
        {
            double[,] matrix = new double[matrixA.GetLength(0), matrixA.GetLength(1)];

            for (int i = 0; i < matrixA.GetLength(0); i++)
            {

                for (int j = 0; j < matrixA.GetLength(1); j++)
                {
                    matrix[i, j] = matrixA[i, j] * konstanta;
                }
            }
            return matrix;
        }
        public double[,] MulOFMatrix(double[,] matrixA, double[,] matrixB)
        {
            double[,] matrix = new double[matrixA.GetLength(0), matrixB.GetLength(1)];

            for (int k = 0; k < matrixB.GetLength(1); k++)
            {
                for (int i = 0; i < matrixA.GetLength(0); i++)
                {
                    for (int j = 0; j < matrixA.GetLength(1); j++)
                    {
                        matrix[i, k] += matrixA[i, j] * matrixB[j, k];
                    }

                }
            }
            return matrix;
        }
        public double[,] Transposition(double[,] matrixA)
        {
            double[,] matrix = new double[matrixA.GetLength(1), matrixA.GetLength(0)];
                for (int i = 0; i < matrixA.GetLength(1); i++)
                {    
                    for (int j = 0; j < matrixA.GetLength(0); j++)
                    {
                        matrix[i, j] = matrixA[j, i];
                    }
            }
            return matrix;
        }
        public double Determinant(double[,] matrix)
        {
            double determinant = 1;
            double[,] det = Copy(matrix);
          det =  GaussElimination(det).Item1.Last();
            for (int i = 0; i < matrix.GetLength(0); i++)
                determinant *= det[i, i];
            return determinant * scalarDet;
        }      
        public double[,] InversionOfMatrix(double[,] matrix)
        {
            if (matrix.GetLength(0) != matrix.GetLength(1))
            {
                Console.WriteLine("matrix neni ctvercova, inverzni matici nelze vypocitat");
                return new double[0, 0];
            }
            double[,] Inversion = new double[matrix.GetLength(0), matrix.GetLength(1)];
            double[,] SplitMatrix = new double[matrix.GetLength(0), matrix.GetLength(1)];
            double[,] MegaMatrix = new double[matrix.GetLength(0), matrix.GetLength(1)*2];
            for (int i = 0; i < Inversion.GetLength(0); i++)
            {
                Inversion[i, i] = 1;
            }
            for (int i = 0; i < matrix.GetLength(0); i++)
                for (int j = 0; j < matrix.GetLength(1); j++)
                    MegaMatrix[i, j] = matrix[i, j];
            for (int i = 0; i < matrix.GetLength(0); i++)
                for (int j = matrix.GetLength(1); j < MegaMatrix.GetLength(1); j++)
                    MegaMatrix[i, j] = Inversion[i, j - matrix.GetLength(1)];
            //spojeni matic
            MegaMatrix = GaussEliminationSupport(MegaMatrix).Item1;
            //rozdeleni matic
            for (int i = 0; i < matrix.GetLength(0); i++)
                for (int j = 0; j < matrix.GetLength(1); j++)
                    SplitMatrix[i, j] = MegaMatrix[i, j];

            for (int i = 0; i < Inversion.GetLength(0); i++)
                for (int j = Inversion.GetLength(1); j < MegaMatrix.GetLength(1); j++)
                    Inversion[i, j - matrix.GetLength(1)] = MegaMatrix[i, j];

            //preklopeni matrix
            double[,] temp1 = new double[SplitMatrix.GetLength(0), SplitMatrix.GetLength(1)];
            double[,] temp2 = new double[Inversion.GetLength(0), Inversion.GetLength(1)];
            for (int i = 0; i < SplitMatrix.GetLength(0); i++)
                for (int j = 0; j < SplitMatrix.GetLength(1); j++)
                    temp1[i, j] = SplitMatrix[SplitMatrix.GetLength(0)-1-i,SplitMatrix.GetLength(0)-1-j];

            for (int i = 0; i < Inversion.GetLength(0); i++)
                for (int j = 0; j < Inversion.GetLength(1); j++)
                    temp2[i, j] = Inversion[Inversion.GetLength(0) - 1 - i, Inversion.GetLength(0) - 1 - j];
            SplitMatrix = Copy(temp1);
            Inversion = Copy(temp2);

            //SpojeniMatic
            for (int i = 0; i < SplitMatrix.GetLength(0); i++)
                for (int j = 0; j < SplitMatrix.GetLength(1); j++)
                    MegaMatrix[i, j] = SplitMatrix[i, j];

            for (int i = 0; i < Inversion.GetLength(0); i++)
                for (int j = matrix.GetLength(1); j < MegaMatrix.GetLength(1); j++)
                    MegaMatrix[i, j] = Inversion[i, j - matrix.GetLength(1)];

            MegaMatrix = GaussEliminationSupport(MegaMatrix).Item1;
            for (int i = 0; i < MegaMatrix.GetLength(0); i++)
            {
                double nasobek = 1/MegaMatrix[i, i];
                for (int j = 0; j < MegaMatrix.GetLength(1); j++)
                    MegaMatrix[i, j] = nasobek*MegaMatrix[i,j];
            }

            for (int i = 0; i < Inversion.GetLength(0); i++)
                for (int j = Inversion.GetLength(1); j < MegaMatrix.GetLength(1); j++)
                    Inversion[i, j - matrix.GetLength(1)] = MegaMatrix[i, j];
           
          
            for (int i = 0; i < Inversion.GetLength(0); i++)
                for (int j = 0; j < Inversion.GetLength(1); j++)
                    temp2[i, j] = Inversion[Inversion.GetLength(0) - 1 - i, Inversion.GetLength(0) - 1 - j];
            Inversion = Copy(temp2);
            double[,] kontrola = MulOFMatrix(matrix, Inversion);
            return Inversion;
        }
        public bool sharpRowDiagonallyDominant(double[,] matrix)
        {
            //podminka konvergence metod
            for (int i = 0; i < matrix.GetLength(0); i++)
            {
                double Sum =0;
                for (int j = 0; j < matrix.GetLength(1); j++)
                {
                    if (i != j)
                    Sum += Math.Abs(matrix[i, j]);

                }
                if (matrix[i, i] < Sum)
                    return false;

            }
            return true;
        }
        public bool sharpColumnDiagonallyDominant(double[,] matrix)
        {
            for (int i = 0; i < matrix.GetLength(0); i++)
            {
                double Sum = 0;
                for (int j = 0; j < matrix.GetLength(1); j++)
                {
                    if (i != j)
                        Sum += Math.Abs(matrix[j,i]);

                }
                if (matrix[i, i] < Sum)
                    return false;

            }
            return true;
        }
        public double EndingConditionOfIteration(double[,] vectorA, double[,] vectorB)
        {
            double[] Array = new double[vectorA.GetLength(0)];
            for (int i = 0; i < vectorA.GetLength(0); i++)
            {
                for (int j = 0; j < vectorA.GetLength(1); j++)
                {
                    Array[i] = Math.Abs(vectorA[i, j] - vectorB[i, j]);
                }
            }
            return Math.Round(Array.Max(), Extension.countDecimalPlaces);
        }

    }
    class LinearniRovnice_ret
    {
        double[,] matrix;
        double[,] matrixL;
        double[,] matrixU;
        double[,] matrixP;
        double[,] vectorB;
        double[,] vectorY;
        double[,] vectorX;
        public double[,] matrixProperty
        {
            get { return matrix; }
            set { matrix = value; }
        }
        public double[,] matrixLProperty
        {
            get { return matrixL; }
            set { matrixL = value; }
        }
        public double[,] matrixUProperty
        {
            get { return matrixU; }
            set { matrixU = value; }
        }
        public double[,] matrixPProperty
        {
            get { return matrixP; }
            set { matrixP = value; }
        }
        public double[,] vectorBProperty
        {
            get { return vectorB; }
            set { vectorB = value; }
        }
        public double[,] vectorYProperty
        {
            get { return vectorY; }
            set { vectorY = value; }
        }
        public double[,] vectorXProperty
        {
            get { return vectorX; }
            set { vectorX = value; }
        }

    }
    class LinearniRovniceSlozky_ret
    {
        List<double[,]> vectors = new List<double[,]>();
        double[,] matrix;
        double[,] vectorB;
        double[,] vectorX;
        List<double> accuracy = new List<double>();
        double factor;
        public double[,] matrixProperty
        {
            get { return matrix; }
            set { matrix = value; }
        }
        public double[,] vectorBProperty
        {
            get { return vectorB; }
            set { vectorB = value; }
        }
        public double[,] vectorXProperty
        {
            get { return vectorX; }
            set { vectorX = value; }
        }
        public List<double[,]> vectorsProperty
        {
            get { return vectors; }
            set { vectors = value; }
        }
        public double factorProperty
        {
            get { return factor; }
            set { factor = value; }
        }
        public List<double> vectorsAccuracyProperty
        {
            get { return accuracy; }
            set { accuracy = value; }
        }
    }

}
